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I. INTRODUCTION 

Monte Carlo methods [H appeared about sixty years 
ago with the need to evaluate numerical values for vari- 
ous complex problems. These methods evolved and were 
applied early to quantum problems, thus putting within 
reach exact numerical solutions to non-trivial quantum 
problems 0, H, 0) [E] ■ Many improvements of these meth- 
ods followed, avoiding critical slowing down near phase 
transitions and allowing to work directly in the continu- 
ous imaginary time limit % 0, i, i, [ul, [H, d . In re- 
cent years, interest in methods that work in the canonical 
ensemble with global updates yet allow access to Green 
functions has intensified [ij, [l^. However, a method 
that works well for a given Hamiltonian often needs ma- 
jor modifications for another. For example, the addition 
of a 4-site ring exchange term in the bosonic Hubbard 
model required special developments for a treatment by 
the stochastic series expansion algorithm flGj, as well as 
by the wordline algorithm [17]. This can result in long 
delays. It is, therefore, advantageous to have at one's 
disposal an algorithm that can be applied to a very wide 
class of Hamiltonians without requiring any changes. In 
a recent publication [l^, the stochastic Green function 
(SGF) algorithm was presented, which meets this goal. 
The algorithm can be applied to any lattice Hamiltonian 
of the form 



n = v-T, 



(1) 



where V is diagonal in the chosen occupation number 
basis and T has only positive matrix elements. This in- 
cludes all kinds of systems that can be treated by other 
methods presented in ref . d, [ll|, [H, [H, [13] , for instance 
Bose-Hubbard models with or without a trap, Bose- Fermi 
mixtures in one dimension, Heisenberg models... In par- 
ticular Hamiltonians for which the non-diagonal part 
T is non-trivial (the eigen-basis is unknown) are easily 
treated, such as the Bose-Hubbard model with ring ex- 
change [Ts, 17], or multi-species Hamiltonians in which a 
given species can be turned into another one (see eq. ([i^]) 
and fig. [3] and |4] for a concrete example). Systems for 
which it is not possible to find a basis in which V is 
diagonal and T has only positive matrix elements are 
said to have a " sign problem" , which usually arises with 



fermionic and frustrated systems. As other QMC meth- 
ods, the SGF algorithm does not solve this problem. 

The algorithm allows to measure several quantities of 
interest, such as the energy, the local density, local com- 
pressibility, density-density correlation functions... In 
particular the winding is sampled and gives access to the 
superfluid density. Equal-time n-body Green functions 
are probably the most interesting quantities that can be 
measured by the algorithm, by giving access to momen- 
tum distribution functions which allow direct compar- 
isons with experiments. All details on measurements are 
given in ref.[l5|. 

In addition the algorithm has the property of being 
easy to code, due in part to a simple update scheme in 
which all moves are accepted with a probability of 1. 
Despite of such generality and simplicity, the algorithm 
might suffer from a reduced efficiency, compared to other 
algorithms in situations where they can be applied. 

The purpose of this paper is to present a "directed" 
update scheme that (i) keeps the simplicity and general- 
ity of the original SGF algorithm, and (ii) enhances its 
efficiency by improving the sampling over the imaginary 
time axis. While the SGF algorithm is not intended to 
compete with the speed of other algorithms, the improv- 
ment resulting from the directed update scheme is re- 
markable (see section V). But what makes the strength 
of the SGF method is that it allows to simulate Hamil- 
tonians that cannot be treated by other methods or that 
would require special developments (see eq. (|49p for a con- 
crete example). The paper is organized as follows: We in- 
troduce in section II the notations and definitions used in 
ref. [l^ . In section III, we propose a simplification of the 
update scheme used in the original SGF algorithm, and 
determine how to satisfy detailed balance. A generaliza- 
tion of the simplified update scheme is presented in sec- 
tion IV, which constitutes the directed updated scheme. 
Finally section V shows how to determine the introduced 
optimization parameters, and presents some tests of the 
algorithm and a comparison with the original version. 



II. DEFINITIONS AND NOTATIONS 

In this section, we recall the expression of the " Green 
operator" introduced in the SGF algorithm, and the ex- 
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tended partition function which is considered. Although 
not required for understanding this paper, we refer the 
reader to ref.[lB] for full details on the algorithm. As 
many QMC algorithms, the SGF algorithm samples the 
partition function 



Z{(3) = Tr e 



(2) 



The algorithm has the property of working in the canon- 
ical ensemble. In order to define the Green operator, we 
first define the "normalized" creation and annihilation 
operators, 



1 



A = 



1 



(3) 



where and a are the usual creation and annihilation 
operators of bosons, and n = a^a is the number operator. 
From ([3]) one can show the following relations for any 
state \nj in the occupation number representation. 



A^\n) = |n + 1) A\n) 



1)^ 



(4) 



with the particular case -4|0) — 0. Appart from this 

exception, the operators and A change a state \rij by 
respectively creating and annihilating one particle, but 
they do not change the norm of the state. 

Using the notation {iplig} to denote two subsets of site 
indices Zi,i2, • • • ,?p and ji,j2, ■ ■ ■ ,jq with the constraint 
that all indices in subset i are different from the indices in 
subset j (but several indices in one subset may be equal) , 
we define the Green operator Q by 



+00 +00 p q 

p=0 q=0 / - I . \ fc=l 1=1 



(5) 



where gpq is a matrix that depends on the application 
of the algorithm [l^. In order to sample the partition 
function an extended partition function Z{(3,t) is 

considered by breaking up the propagator e"''^, and in- 
troducing the Green operator between the broken parts, 

r) = Tr e-^'^-^^^^e"^^. (6) 

Defining the time dependant operators T(t) and G{t), 

f{T) = e-^^fe-^^ g{T) e^^Ge-^^, (7) 

and working in the occupation number basis in which V is 
diagonal, the extended partition function takes the form 

Z(/3,t)=^ /(V'o|e-''^r(T„)|V'„_i)(V'„-i|r(r„_i)|V'„_2) 

„>oio<ri<-<r„</3 

X • • • (V'L+i |r(Ti) I I V'fl,)(V'i?Jt(Tfl.)| V-fl-i) (8) 

X • • • {lp2\T{T2)\lpl){lJjl\'T{Ti)\xpQ)dTi ■ ■ -dTn, 



where the sum J2n>o implicitly runs over complete sets 
of states {jV'fc)}- We will systematically use the labels 
L and R to denote the states appearing on the left and 
the right of the Green operator, and use the notation T4 
to denote the diagonal energy (V'fcj VjV'fc)- We will also 
denote by tl and tr the time indices of the T operators 
appearing on the left and the right of G- 

As a result, the extended partition function is a sum 
over all possible configurations, each being determined 
by a set of time indices ti, T2, • • • , t/j, t, t^, • • • , t„ and a 
set of states |V'o), ■ ■ ■\^r),\'4'l), ■■■\i^n~i)- The 

algorithm consists in updating those configurations by 
making use of the Green operator. Assuming that the 
Green operator is acting at time r, it can "create" a T 
operator (that is to say a T operator can be inserted in 
the operator string) at the same time, thus introducing 
a new intermediate state, then it can be shifted to a dif- 
ferent time. While shifting, any T operator encountered 
by the Green operator is "destroyed" (that is to say re- 
moved from the operator string). Assuming a left (or 
right) move, creating an operator will update the state 
V'j?) (or \iPl)), while destroying will update the state 
iPl) (or I '(/'-«))■ When a diagonal configuration of the 
Green operator occurs, = [V'i?); such a configura- 

tion associated to the extended partition function ([5]) is 
also a configuration associated to the partition function 
(12). Measurements can be done when this occurs (see 
ref. [15!] for details on measurements) . 

Next section presents a simple update scheme that 
meets the requirements of ergodicity and detailed bal- 
ance. 



III. SIMPLIFIED UPDATE SCHEME 

Before introducing the directed update, we start by 
simplifying the update scheme used in the original SGF 
algorithm. 



A. The update scheme 

We will assume in the following that a left move of 
the Green operator is chosen. In the original version, the 
Green operator G{t) can choose to create or not on its 
right a T operator at time r. Then a time shift At to 
the left is chosen for the Green operator with an expo- 
nential distribution in the range [0; -|-oo[. If an operator 
is encountered while shifting the Green operator, then 
the operator is destroyed and the move stops there. As 
a result, four possible situations can occur during one 
move: 



1. No creation, shift, no destruction. 

2. Creation, shift, no destruction. 

3. No creation, shift, destruction. 
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4. Creation, shift, destruction. 

It appears that the first possibihty "no creation, no de- 
struction" is actually useless, since no change is per- 
formed in the operator string. The idea is to get rid of 
this possibility by forcing the Green operator to destroy 
an operator if no creation is chosen. A further simplifi- 
cation can be done by noticing that the last possibility 
"creation, destruction" is not necessary for the ergodicity 
of the algorithm, and can be avoided by restricting the 
range of the time shift after having created an operator. 
Therefore we replace the original update scheme by the 
following: We assume that the Green operator is acting 
at time r and that the operator on its left is acting at 
time TL- The Green operator Q{t) chooses to create or 
not an operator on its right at time r. If creation is cho- 
sen, then a time shift At of the Green operator is chosen 
to the left in the range [0; tl — r[, with the probability 
distribution defined below. If no creation is chosen, then 
the Green operator is directly shifted to the operator on 
its left at time r^, and the operator is destroyed. As a 
result only two possibilities have to be considered: 

1. Creation, shift. 

2. Shift, destruction. 

Figure [1] shows the associated organigram. Section III.B 
explains how detailed balance can be satisfied with this 
simplified update scheme. 




V 

End 



FIG. 1: The simplified update scheme. See text for details. 



B. Detailed balance 

When updating the configurations according to the 
chosen update scheme, we need to generate different tran- 
sitions from initial to final states with probabilities that 
satisfy detailed balance. In this section we propose a 
choice for these probabilities, and determine the corre- 
sponding acceptance factors. We denote the probability 
of the initial (final) configuration by Pi (Pf). We denote 
by Si^f the probability of the transition from configu- 
ration i to configuration /, and by S the probability 
of the reverse transition. Finally we denote by Ai^f the 
acceptance rate of the transition from i to f, and by 
A the acceptance rate of the reverse transition. The 
detailed balance can be written as 

PA^fA^f - PfSf^.Af^,. (9) 

We will make use of the Metropolis solution [Tsj . 

^,^/ =min(l,<7) (10) 

with 



We will use primed (non-primed) labels for states and 
time indices to denote final (initial) configurations. 

1. Creation, shift 

We consider here the case where a left move is chosen, 
an operator is created on the right of the Green operator 
at time r, and a new state is chosen. Then a time shift 
to the left is chosen for the Green operator in the range 
[0, — r}j[. It is important to note that and corre- 
spond to the time indices of the operators appearing on 
the left and the right of the Green operator after the new 
operator has been inserted, that is to say at the moment 
where the time shift needs to be performed. Thus we 
have = tl and t'j^ — t. 

The probability of the initial configuration is the Boltz- 
mann weight appearing in the extended partition func- 
tion dS]): 

P^ CX {ijL\G{T)\ijR) 

oc e^^^(VL|^|^i?>-"^^ (12) 
The probability of the final configuration takes the form: 
Pf CX (^l|^(r')|V'k)('Ak|^(T;?)|V'k_i) 

It is important here to realize that the Green operator 
only inserted on its right the operator I^Pr) {i'nl'i' {T'ft) , 
before being shifted from to r'. Therefore we have the 
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equalities (V^lI = (V'l], = IV'fl), "^l = ^l, and 

The probability Si^f of the transition from the initial 
configuration to the final configuration is the probability 
P{^) of a left move, times the probability -PJ_(t) of a 
creation, times the probability P^{ip'fi) to choose the new 
state V'fl, times the probability i'^' ~'^'b) to shift the 
Green operator by t' — t'j^, knowing that the states on the 
left and the right of the Green operator at the moment 
of the shift are ■0^ and ip'j^: 

S,^j = P{^)pI{t)P^{^'^)pE"' {r' - t'^) (14) 

The probability of the reverse transition is simply the 
probability P{-^') of a right move, times the probability 
of no creation, 1 — P1^{t'): 



Sf^. = P{-^')[l-Pl{r')] 



(15) 



From the original version of the SGF algorithm, we know 
that choosing the time shift with an exponential distribu- 
tion is a good choice, because it cancels the exponentials 
appearing in the probabilities of the initial (|12p and final 
(I13p configurations, avoiding exponentially small accep- 
tance factors. However a different normalization must 
be used here, since the time shift is chosen in the range 
[0;t^ — t'jj[ instead of [0; +oo[. The suitable solution is: 



P£^'(Ar) 



1 -e-(^i-^«)(^A-^^i) 



(16) 



It is straightforward to check that the above probabil- 
ity is correctly normalized and well-defined for any real 
value of V^ — Vl, the particular case VI — reducing to 
the uniform distribution P{At) = 1/ (r^ — r^) (note that 
— t|j is always a positive number). For the probabil- 
ity P^{tlj'j^) to choose the new state ip'j^, the convenient 
solution is the same as in the original version: 



{iPl\GT\-iPr) 



(17) 



Putting everything together, the acceptance factor (fTTj) 
becomes 





\5f\i'R) 






^)Pl(r) 



Pi-^')[1-Plir')] [l- 



where we have used the notation to emphasize that 
this acceptance factor corresponds to a creation. It is 
also important for the remaining of this paper to note 
that q1_ is written as a quantity that depends on the 
initial configuration, times a quantity that depends on 
the final configuration. 



2. Shift, destruction 

We consider here the case where a left move is cho- 
sen, and the operator on the left of the Green operator is 
destroyed. This move corresponds to the inverse of the 
above "creation, shift" move. Thus, the corresponding 
acceptance factor gf, is obtained by inverting the accep- 
tance factor g^, exchanging the initial time r and final 
time r', and switching the direction. However tl — tj^ 
represents an absolute time shift, so tl and th do not 
have to be exchanged. We get 

d ^ Vl-Vr 

P(^) [1 - pI{t)] [1 - e-(^^-^«)(^^-^«)] 



{^'lI-PSI^'r) 



(19) 



which is written as a quantity that depends on the initial 
configuration, times a quantity that depends on the final 
configuration. 



3. Simplification of the acceptance factors 

We will use here the short notation (^), {GP), and 

(PG) to denote respectively the quantities (fpLlGlipR), 

(VlI^^IV't?); and (iPlIPGIiPr)- As in ref. [l^, we have 
some freedom for the choice of the probabilities of choos- 
ing a left or right move, P(^) and P(^), and the prob- 
abilities of creation PJL(t) and P1^{t). A suitable choice 
for those probabilities can be done in order to accept all 
moves, resulting in an appreciable simplification of the 
algorithm. For this purpose, we impose the acceptance 
factor (or q'i^) to be equal to the acceptance factor 
gf_ (or q'^). This allows to determine the probabilities 
PI{t) and P1(t), 



p1(t) = 
Plir) = 



(GP) 



(PG) + (G) Y—^iTZ^T^^v^rrvJJ 



(20) 



, (21) 



and the acceptance factors = and q't^ = q*^ take 
the form 



P(^)r^(r') 



P{^)rAr'y 



with 



Vl~Vr 



(St) 

^ 1 - e-(ri:-rH)(yi,-VR) 

Vr~Vl 



(23) 
(24) 
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Finally we can impose the acceptance factors and 
to be equal. This implies 



r^(T) + r^(T) 

(25) 

Defining i?(T) = r^{T) + r^(r), we are left with a single 
acceptance factor, 



Rjr) 



(26) 



which is independent of the chosen direction, and inde- 
pendent of the nature of the move (creation or destruc- 
tion). Thus all moves can be accepted by making use 
of a proper reweighting, as explained in ref. p^ . The 
appendix shows how to generate random numbers with 
the appropriate exponential distribution ((T| 



C. Discussion 

Although the above simplified update scheme works, 
it turns out to have a poor efficiency. This is because 
of a lack of "directionality": The Green operator has, 
in average, a probability of 1/2 to choose a left move or 
a right move. Therefore the Green operator propagates 
along the operator string like a " drunk man" , with a 
diffusion-like law. The basic creation and destruction 
processes correspond to the steps of the random walk. 

This suggests that the efficiency of the update scheme 
can be improved if one can force the Green operator to 
move in the same direction for several iterations. Next 
section presents a modified version of the simplified up- 
date scheme, which allows to control the mean length of 
the steps of the random walk, that is to say the mean 
number of creations and destructions in a given direc- 
tion. The proposed directed update scheme can be con- 
sidered analogous to the "directed loop u pda te" used in 
the stochastic series expansion algorithm jlll . [T9| , which 
prevents a worm from going backwards. However the 
connection should not be pushed too far. Indeed the pic- 
ture of a worm whose head is evolving both in space and 
imaginary time accross vertices is obvious in a loop algo- 
rithm. In such algorithm, a creation (or an annihilation) 
operator which is represented by the head of a worm is 
propagated both in space and imaginary time, while an 
annihilation (or a creation) operator represented by the 
tail of the worm remains at rest. The loop ends when the 
head of the worm bites the tail. 

Such a worm picture is not obvious in the SGF algo- 
rithm: Instead of single creation or annihilation opera- 
tors, it is the full Green operator over the whole space 
that is propagated only in imaginary time. This creates 
open worldlines, thus introducing discontinuities. These 
discontinuities increase or decrease while propagating in 
imaginary time. All open ends of the worldlines are lo- 
calized at the same imaginary time index. Therefore it is 
actually not possible to draw step by step a worm whose 



head is evolving in space and imaginary time until it bites 
its tail. 



IV. DIRECTED UPDATE SCHEME 

We present in this section a directed update scheme 
which is obtained by modifying slightly the simplified up- 
date scheme, thus keeping the simplicity and generality 
of the algorithm. 



A. The update scheme 

Assuming that a left move is chosen, the Green opera- 
tor chooses between starting the move by a creation or a 
destruction. After having created (or destroyed) an op- 
erator, the Green operator can choose to keep moving in 
the same direction and destroy (or create) with a proba- 
bility (or P^L^), or to stop. If it keeps moving, then a 
destruction (or creation) occurs, and the Green operator 
can choose to keep moving and create (or destroy) with 
a probability (or P^)... and so on, until it decides 
to stop. If the last action of the move is a creation, then 
a time shift is chosen. The organigram is represented in 
Figure H 



B. Detailed balance 

In order to satisfy detailed balance, in addition to 
the acceptance factors q'l_ and gf,, we need to deter- 
mine new acceptance factors of the form g^cdcdc -- 
qdcdcdcd- -^ We first determine the new expressions of q1_ 
and gf_ resulting from the directed update scheme. For 
g^, the previous probability Si^f has to be multiplied 
by the probability to stop the move after having created, 
1 — P^{t'). The previous probability Sf^i has to be 
multiplied by the probability to stop the move after hav- 
ing destroyed, 1 — P'^ (r) . We get for q1_ and q"^ the new 
expressions: 



at 



(^L|gr|^fl)[i-p^'=(T)] 

{M5\^PR)Pi^)PliT) 

P(^')[l - Pli^')] [1 - e-(^^.-^«)(^H-^^)" 
[1-P_^'^(r)] (Vl-Vr) 



■27) 



P(^)[l - PI{t)] [l - e-(^^-^«)(^^-^«)] 
{ij'^\fg\^'^)[l-PiL^{T')y 



(28) 



1. Creation, destruction 

We consider here the case where a left move is chosen, 
an operator is created on the right of the Green opera- 
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FIG. 2: The directed update scheme. See text for details. 

tor, and a new state is chosen. Then the operator on the 
left of the Green operator is destroyed. Using the super- 
scripts a, 6, c, • • • to denote intermediate configurations 
between initial and final configurations, the sequence is 
the following 

1. Pi OC (V'L+l|t(TL)|VL)<V'i|^W|V'fl) 

2. (V'£+i|7-(T£)|V'£)(V£|^(T")|V'^)(V'^|r(T]|)|V'S-i) 

where we have = {ipt+i] — (V'lIj 

i'R) = \rR-i) = |V'l)(V'l| = and 

''Pr){'^r\ ~ • The probability of the transi- 

tion from the initial configuration to the final configura- 
tion is the probability f (•«— ) to choose a left move, times 
the probability PJL(t) to create an operator at time t, 
times the probability P^{ip%) to choose the new state 
V'^, times the probability P^{a) to keep moving and de- 
stroy, times the probability 1 — P^{t') to stop the move 
after having destroyed: 

S^^^f = P{^)Pl{T)P^{rR)P^J{a)[l - Pt^ir')] (29) 

The probability of the reverse move is exactly symmetric: 

Sf^i = P{^')Pl{T')P^{rL)Pt'{a) [1 - Pt%T)] (30) 



It is important to notice that, when in the intermediate 
configuration a, the time r£ of the operator to the left of 
the Green operator is equal to tl , and the time of the 
operator to the right of the Green operator is equal to r. 
Thus the acceptance factor takes the form 

, ^ {,PL\gf\i'R)[l-Pt%r)] 

{^PL\gMp{^)pl{T) 

__1 L 

e-(^^-^s)^i"pfe<i(a) 
^ {^y,\g\^'^)P{-^')Pl{r') 

{i,i\fg\i,'^)[i-P!^^ir')y 

and is written as a quantity that depends on the initial 
configuration, times a quantity that depends on the inter- 
mediate configuration a, times a quantity that depends 
on the final configuration. It is useful for the remaining 
of the paper to define the intermediate acceptance factor, 

q^-'^ia) = ^ (32) 

2. Destruction, creation 

We consider here the case where a left move is cho- 
sen, the operator on the left of the Green operator is 
destroyed, then an operator is created on its right, and a 
new state is chosen. Finally a time shift is chosen. The 
sequence of configurations is the following 

1. P, CX {i'L+l\fiTL)\^L){^L\g{T)\^R) 

2. (V2|^(r'')|V^) 

3. Pf OC (Vi|^(T')|V;j)(V'k|^(^ii)|V';j-i), 

where we have = — (V't|) and 

IV'-r) = IV'fl) = li^R-i)- The probability of the tran- 
sition from the initial configuration to the final config- 
uration is the probability P{^) to choose a left move, 
times the probability 1 — Pl_ (t) of no creation, times the 
probability P^{a) to keep moving and create, times the 
probability P^{ij/j^) to choose the new state ij/j^, times 
the probability 1 — P^{t') to stop the move after having 
destroyed, times the probability ^ (r' — rjj) to shift 
the Green operator by r' — r^: 

Si^f = Pi^)[l-Pl{T)]Pt%a)P^{^'^) 

X [1-P^{t')]pE^'{t' -t'h) (33) 

The probability of the reverse move is exactly symmetric: 

Sf^i = Pi-^')[l-Plir')]Pt%a)P^{^i^) 

X [l-P^f(T)]Pi:^(TL-T) (34) 
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The acceptance factor takes the form 



[1 - Pi'^(T)] [Vl - Vr) 



P(^) [1 - Pt(T)] [1 - e-(^^--^«)(^^--^«)] 

{rL\G'T\rR)pt'{a) 



{'il;l\fg\rR)Pt^{a) 
P{^')[1-PI{t')\ [1 



r'H)(vk-vL)^ 



[1 _ p^d^r')] (y^ 



and is written as a quantity that depends on the initial 
configuration, times a quantity that depends on the inter- 
mediate configuration a. times a quantity that depends 
on the final configuration. It is useful for the remaining 
of the paper to define the intermediate acceptance factor, 



\a) = 



(^£|gf|V;^)P^^(a) 
{rL\TQ\rR)Pt{a)' 



(36) 



4. Pf cx (Vi|a(T')|^k)<Vk|^(^fl)|^k-i), 

Considering the intermediate configurations a and b be- 
tween the intial and final configurations, it is easy to show 
that the corresponding acceptance factor can be written 



Generalization 



(38) 



It is straighforward to show that the acceptance fac- 
tors of the form g^^dc^ ^cdcdcdc^ qcdcdcdcdc.. . ^dcdcd^ 

qdcdcdcd^ ^dcdcdcdcd , , .-j CXprCSScd aS products of 

the acceptance factor q1_ (or q'l_) and the intermediate 
factors q^'^ and q'^'^. 

In the same manner, the acceptance factors of the form 

^cdcd ^cdcdcd ^cdcdcdcd ^ , , ^dcdc ^dcdcdc ^dcdcdcdc ^ , \ 

y* — 1 y* — 7 Hi — * * * y* — ? y* — ? y* — * * */ 
can be expressed as products of the acceptance factor 
(or q"^) and the intermediate factors q'^'^ and gfr°. 



3. Creation, destruction, creation 



6. Simplification of the acceptance factors 



We consider here the case where a left move is chosen, 
an operator is created on the right of the Green operator, 
then the operator on its left is destroyed, then a second 
operator is created on its right. Finally, a time shift 
of the Green operator is performed. The sequence of 
configurations is the following 

1. Pi cx <VL+l|t(ri)|VL)(VL|^(T)|V'fl) 

2. (V'£+i|t(TE)|V£)(V'!|a(T«)|^^)(^^|r(r]|)|^^_i) 

3. mQ{r')mmr{ri)\i^'R.,^ 

4. Pf <x (Vi|^(r')|^k)(Vk|^(^K)kk-i)(^k-i|^(T;i-i) 

Considering the intermediate configurations a and h be- 
tween the intial and final configurations, it is easy to show 
that the corresponding acceptance factor can be written 



Here again it is possible to take advantage of the free- 
dom that we have for the choice of the probabilities 
P(^), Pi, P^^ and P^ (or P(^), Pi,, P^ , and P'^). 
A proper choice of these probabilities can be done in 
order to allow us to accept all moves, simplicity and gen- 
erality being the leitmotiv of the SGF algorithm. 

For this purpose, we impose to all acceptance factors 
corresponding to left (or right) moves to be equal. This 
requires the intermediate acceptance factors q'^'^ and 
qd-c i^yj. qc-d gf7°) to be equal to 1. This is re- 
alized if 



pkc _ 



pkc _ 



^cdc 



qt X q^^^ia) x qt%b) 



(37) 



pkd 
pkd 



a.mm,l,^j 

aemm,l,^j 

a,min(l,e-(^'-^«)K'-^-')) 
arfmin(l,e"(^^S-^«)(^^-^«)), 



(39) 

(40) 

(41) 
(42) 



^. Destruction, creation, destruction 



where ac and ad are optimization parameters belonging 
to [O; 1[. By tuning these parameters, the mean length of 
the steps of the Green operator can be controlled. Note 
that we have explicitly excluded 1 from the allowed val- 
ues for these optimization parameters. This is necessary 
for the Green operator to have a chance to end in a di- 
agonal configuration, I'tpL) = |'0-r)- Indeed, the choice 
etc = cud = i would systematically lead to values of 1 for 
1. Pi oc <VL+2|t(rL+i)|^L+i)(^L+i|t(Ti)|^i)(^i|a(r)feprobabilitiesP'=^andP'='^ for diagonal configurations. 

Therefore the Green operator would never stop in a diag- 
onal configution, and no measurement could be done. It 
is important here to note that the quantities {G}, {G'T}, 
and {'TQ) are evaluated between the states on the left 



We consider here the case where a left move is cho- 
sen, the operator on the left of the Green operator is 
destroyed, then an operator is created on its right. Fi- 
nally a second operator on the left of Green operator is 
destroyed. The sequence of configurations is the follow- 
ing 



2. {rL+i\f{Ti)\rL){rL\GiT^)\rR) 

3- (V'i+i|t(ri)|V'i)(^i|^(r^)|V^)(^|^|r(ri)|V'|,_,) 
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and the right of the Green operator that are present at 
the moment where those quantities are needed, as well 
as for the times indices and and the potentials Vj^ 
and V^. 

All acceptance factors corresponding to a given direc- 
tion of propagation become equal if we choose for the 
creation probabilities: 



Plir) 



Plir) 



{Vl-Vr) 



-443) 



(T-G) 



{fG) + {G) 



{Vr~Vl) 



Finally, all acceptances factors become indepen- 
dant of the direction of propagation if we choose 

— i^xl'^ I \ with 



P(-) = and PH) 



= [i-Pi1 



r^{T) = [I - Pt']^^^^ 



[1 - Pt^] {Vl ~ Vr) 
[l-Pt^]iVB.-VL) 



-(45) 



.(46) 



As a result all moves can be accepted again, ensuring the 
maximum of simplicity of the algorithm. We still have 
some freedom for the choice of the optimization parame- 
ters ac and aj,. This is discussed in next section. 



seems reasonable to impose a condition of uniform sam- 
pling, P'^'^ = p^'^ , This condition can be satisfied by ad- 
justing dynamically the values of etc and ad during the 
thermalization process. For this purpose we introduce 
a new optimization parameter a € [O; 1 [ and apply the 
following algorithm from time to time while thermalizing 
(we start with ac — Ud — oi): 

Evaluate P''" and P^"^ over few iterations 

j£ pkc ^ pkd 

pkc 

'■pkd 
pkd 



then ad ^ ad- 



else ac — 
If ttc < ad 



'^'^ pkc 



then ac = — , ad 
ad 

else ad = — , ac = 
ac 



Thus we are left with the optimization parameter a. In 
order to determine the optimal value, we have considered 
2 different Hamiltonians TLi and 7^2, and evaluated the 
efficiency of the algorithm while scanning a. The first 
Hamiltonian we have considered describes free hardcore 
bosons and is exactly solvable. 



Hi 



-t ( a\aj+a]a. 



(48) 



V. TEST AND OPTIMIZATION OF THE 
ALGORITHM 

From the central limit theorem, we know that the er- 
rorbar associated to any measured quantity must de- 
crease as the square root of the number of measurements, 
or equivalently, the square root of the time of the simu- 
lation. Therefore it makes sense to define the efficiency 
£^ of a QMC algorithm by 



£{n,o) = 



1 



T(f2) Ae'(n) 



(47) 



where VL represents the set of all optimization parameters 
of the algorithm, O is the measured quantity of interest, 
T{Vt) is the time of the simulation, and AO (SI) is the 
errorbar associated to the measured quantity O. This 
definition ensures that £ is independent of the time of the 
simulation. As a result, the larger £ the more efficient 
the algorithm. 

In the present case we have f2 = |ac, a^}, while Q, — % 
for the original SGF algorithm. It is useful here to real- 
ize that, by symmetry, the mean values of P^ and 

(and P^ and P^'*) must be equal. Therefore we define 
pkc ^ (^pkc\^ ^ ^pfec^ g^^^ pkd ^ ^pkd\^ ^ {P^). It 



where the sum runs over pairs of first neighboring sites 
and t is the hopping parameter. The second Hamiltonian 
is highly non-trivial and describes a mixture of atoms and 
diatomic molecules, with a special term allowing conver- 
sions between the two species [20j . 



^2 = 



Uaa n - (n? - 1) + Umm ^ ("f " l) + Ua„^ 
i i 



where a| and (mj and m^) are the creation and an- 
nihilation operators of atoms (molecules), t^, tm, Uaa, 
Umm, and Uam are respectively the hopping parameter 
of atoms, the hopping parameter of molecules, the atomic 
onsite interaction parameter, the molecular onsite inter- 
action parameter, and the inter-species interaction pa- 
rameter. The conversion term is tunable via the param- 
eter g and does not conserve the number Na of atoms or 
the number Nm of molecules. However the total num- 
ber of particles N = Na + 2Nm is conserved and is the 
canonical constraint. The parameter D allows to control 
the ratio between the number of atoms and molecules. 
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a 


(S(a)) £{a,E) 


£{a,ps) 


£{a,n{Q)) 


a 


£{a,na{0)) £[a,n^{Q)) 


£{a,Va) 


£(a,Vm) 





1.00 


0.307400 


0.487457 


0.503105 





0.433382 


0.234412 


1.323720 


0.239113 


0.1 


1.10 


0.774161 


0.513633 


0.805048 


0.1 


0.269700 


0.181019 


0.585183 


0.248060 


0.5 


1.91 


0.430843 


3.771422 


1.289757 


0.5 


1.752466 


2.806166 


2.667114 


1.357462 


0.9 


7.00 


0.977413 


5.400997 


6.629893 


0.9 


7.080124 


5.638859 


16.454676 


4.482435 


0.95 


10.49 


2.427874 


10.688100 


7.994883 


0.95 


4.893878 


3.757436 


5.088775 


2.248427 


0.99 


17.49 


1.286403 


27.281408 


1.327064 


0.99 


3.871723 


2.341222 


7.783268 


1.279447 


0.9999 


20.93 


0.818048 


17.510068 


1.059823 












0.999999 


21.00 


0.710448 


13.353809 


0.779245 













TABLE I: Relative efRciency of the algorithm applied to Tii 
at half filling for the energy, the superfluid density, and the 
number of particles in the zero momentum state. 



TABLE IIL Relative efRciency of the algorithm applied to 
7t!2 for the occupation of the zero momentum state and the 
visibility of atoms and molecules. 



{S{a)) £{a,E) £{a,Pa) £{a,Pn.) 





0.1 
0.5 
0.9 
0.95 
0.99 



1.00 1.086334 0.455569 1.670239 

1.10 1.424308 0.506873 1.858339 

1.88 2.813905 1.265620 4.640123 

6.35 2.562529 5.999027 21.993900 

8.99 2.335315 3.917233 14.361774 

12.79 2.592328 1.721519 6.311612 



TABLE IL Relative efficiency of the algorithm applied to 7^2 
for the energy, and the density of atoms and molecules. 



The application of the SGF algorithm to the Hamilto- 
nian is described in details in ref.[lH|. The changes 
coming with the directed update scheme are completely 
independent of the chosen Hamiltonian. 

The following table shows the mean number of cre- 
ations and destructions in one step, (5'(a)), and the rel- 
ative efRciency f (a, 0)/£{%, O) of the algorithm apphed 
to Til at half filling, for which we have measured the 
energy E, the superfluid density ps, and the number of 
particles in the zero momentum state n{k = 0): 

For 712, we have used the parameters ta = l,tm = 1/2, 
Uaa = 5, Umm = 5, Uam = 5, g = 5, D = 3, and a density 
of particles p = 2. The following tables shows (5'(q!)), 
and the relative efficiency of the algorithm for the energy 
E, the density of atoms and molecules pa and p„i, the 
occupation of the zero momentum state for atoms and 
molecules na{Q) and nm(0), and the atomic and molecu- 
lar visibilities Vq and Vm- 

While the best value of a depends on the Hamiltonian 
which is considered and the measured quantity, it appears 
that a good compromise is to choose a between 0.90 and 
0.99. The improvment of the efficiency is remarkable. 
In the following, we illustrate the applicability of the 
algorithm to problems with non-uniform potentials, by 
adding a parabolic trap to the Hamiltonian 



The parameters Wa and Wm allow to control the cur- 
vature of the trap associated to atoms and molecules, 
respectively, and L is the number of lattice sites. The 
inclusion of this term in the algorithm is trivial since 
only the values of the diagonal energies Vl and Vr are 
changed. Figures ([3]) and (|4]) show the density profiles 
and momentum distribution functions obtained for a sys- 
tem with L = 70 lattice sites initially loaded with 50 
atoms and no molecules, and the parameters ia = 1, 

t„ = 0.5, Uaa = 4, Uam = 12, C/mm = CO , g = 0.5, D = {), 

Wa = 0.008, Wm = 0.008, and /? = 20. The presented 
results have been obtained by performing 10^ updates 
for thermalization, and 2 x 10^ updates with measure- 
ments (an update is to be understood as the occurence 
of a diagonal configuration) . The time of the simulation 
is about 8 hours on a cheap 32 bits laptop with IGHz 
processor, with an implementation of the algorithm in- 
volving dynamical structures with pointers (see ref . [lH ) . 



0,75 



13 0,5 
o 



0,25 
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• • Atoms 

■— ■ Molecules 
' — ' Atoms+Molecules 



60 



70 



nT = WaY,i^-LI2f 



Wr^ 



.Y,{i-Ll2) 



FIG. 3: (Color online) An example of density profiles ob- 
tained when adding the trapping potential (|50p to the Hamil- 
tonian ((49]). The errorbars are smaller than the symbol sizes, 
and are the biggest in the neighborhood of site indices 23 and 
47 where they equal the size of the symbols. 



(50) 
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FIG. 4: (Color online) An example of momentum distribution 
functions obtained when adding the trapping potential (|50[) 
to the Hamiltonian (|49|) . The errorbars are smaller than the 
symbol sizes, and are the biggest for A; = where they equal 
the size of the symbols. 



we have at our disposal a uniform random number gen- 
erator that generates a random variable U with the dis- 
tribution puW) ~ 1 for u E [O; l[, we would like to find 
a function / such that the random variable T = f{U) is 
generated with the distribution 



At.AV, ^ 

Pt ' (V 



AVe- 



-tAV 



,-AtAV 



re[0;Ar[, (51) 



where At and AV^ are the parameters of the exponential 
distribution. Because of the relation T = f{U), the prob- 
ability to find T in the range [t; t -|- dr [ must be equal to 
the probability to find U in the range [u; u + du[. This 
implies the condition 



pu{u)\du\ 



At,AV / \\ 1 I 

Pt [V\dT\, 



(52) 



VI. CONCLUSION 



with 1^1 

I du I 



±#. Thus we have 

du 



We have presented a directed update scheme for the 
SGF algorithm, which has the properties of keeping the 
simplicity and generality of the original algorithm, and 
improves significantly its efficiency. 
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A. Appendix: Exponential random number 
generator 

We describe here how to generate numbers with the ap- 
propriate exponential distribution p6p . Assuming that 



^yg-/(u)Ay 

1 „ e-ArAV ^„ 



±1. 



(53) 



Taking the anti-derivative with respect to u on both sides 
of the equation, we get 



_p-f{u)AV 



1 



-AtAV 



±{u + C), 



(54) 



where C is a constant. This constant and the correct sign 
are determined by imposing the conditions /(O) = and 
/(I) = At. As a result, if u is a realization of U, then a 
realization of T is given by 



AV 



In [l-u(l-e-^"^^)]. 



(55) 
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